Phases and magnetization process of an anisotropic Shastry-Sutherland model 
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We examine ground state properties of the spin- 1/2 easy-axis Heisenberg model on the Shastry- 
Sutherland lattice with ferromagnetic transverse spin exchange using quantum Monte Carlo and 
degenerate perturbation theory. For vanishing transverse exchange, the model reduces to an antifer- 
romagnetic Ising model that besides Neel order harbors regions of extensive ground state degeneracy. 
In the quantum regime, we find a dimerized phase of triplet states, separated from the Neel ordered 
phase by a superfluid. The quantum phase transitions between these phases are characterized. The 
magnetization process shows a magnetization plateau at 1/3 of the saturation value, that persists 
down to the Ising limit, and a further plateau at 1/2 only in the quantum regime. For both plateaus, 
we determine the crystalline patterns of the localized triplet excitations. No further plateaus nor 
supersolid phases are found in this model. 
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I. INTRODUCTION 

Quantum magnets exhibit a wealth of interesting phe- 
nomena, in particular on low-dimensional frustrated lat- 
tices, where both enhanced quantum fluctuations and ge- 
ometric frustration can destroy semi-classical magnetic 
order. Indeed, various compounds have been character- 
ized to provide realizations of the above paradigm. Re- 
cent examples include the valence bond solids found in 
(C 2 H5)(CH 3 ) 3 P[Pd(dmit) 2 ] 2 i and Zn.Cu^^OD^CkA 
Another material, that has been intensively studied is 
SrCuB 2 (B0 3 ) 2 (we refer to Ref. 1 for a detailed review 
of the various experimental and theoretical explorations 
on this system) . This compound is described well by the 
dimer singlet ground state proven exactly previously by 
Shastry and Sutherland to exist4 in the spin- 1/2 Heisen- 
berg model on the orthogonal dimer lattice, shown in Fig. 
1 . Recently, new results on the magnetization process of 
SrCuB 2 (BC>3) 2 have been presented. In particular, mag- 
netization plateaus at 1/5, 1/6, 1/7, 1/9, and 2/9 of the 
saturated magnetization have been reported^, in addi- 
tion to the previously established plateaus at 1/8, 1/4 
and 1/3. While the existence of some of the reported 
plateaus is at the moment controversial, new B 11 NMR 
data — and magnetic torque measurements^ provide evi- 
dence in favor of a persistent crystalline structure of mag- 
netic excitations also above the 1/8 plateau. The pres- 
ence of intra-dimer Dzyaloshinskii-Moriya interactions in 
SrCuB 2 (B03) 2 however calls for a more complex scenario 
than a direct interpretation in terms of supersolidity of 
triplet excitations in this regime^. 

As another realization of the Shastry-Sutherland ge- 
ometry the rare earth tetraborid TmB^ has recently 
been studied in finite magnetic fields. In contrast to 
SrCuB 2 (BC>3) 2 , this metallic compound exhibits stable 
long-range antiferromagnetic order in zero field below 
about 9.8 K. Since full saturation can be obtained for 
magnetic fields parallel to the c-axis above 5T, T111B4 
allows for a complete scan of its magnetization process. 
For this compound magnetization plateaus have been ob- 
served e.g. at fractions 1/2, 1/7, 1/8 and 1/9 of mag- 
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FIG. 1: The orthogonal dimer lattice of the Shastry- 
Sutherland model with spin- 1/2 degrees of freedom on the 
square lattice vertices, and intra-dimer coupling J' (solid 
lines) and nearest neighbor (inter-dimer) coupling J (dashed 
lines). 



netic saturation. Despite the metallic nature of T111B4, 
its magnetism has been suggested to realize an easy-axis 
anisotropic version of the Shastry-Sutherland model close 
to the Ising limit with similar intra- and inter-dimer cou- 
pling strengths 9 . 

In light of the progress in realizing novel quantum 
phases in frustrated quantum magnets, it is important 
to explore in detail the interplay between geometric frus- 
tration and quantum fluctuations in such systems based 
on effective spin models. In many aspects, numerical 
studies have become especially important as an unbiased 
approach to quantum magnetism. However, numerical 
studies of even simple models of frustrated spin systems 
suffer from severe restrictions on the finite sizes accessible 
to current simulation techniques. In particular, quantum 
Monte Carlo (QMC) simulations are tampered by a no- 
torious sign-problemi^ due to odd-length spin-exchange 
paths appearing on non-bipartite lattices. This usually 
restricts unbiased numerical studies to the small lattices 
accessible to exact numerical diagonalization. Notewor- 
thy in this respect are however recent studies employing 
the density matrix renormalization group algorithm on 



2 



the triangular and kagomc lattice Heisenberg mode l 11 ! 12 . 

Here, we employ a different approach in order to ex- 
plore the interplay between quantum fluctuations and 
frustration, by studying a model of quantum magnetism 
in a parameter regime, where geometric frustration is re- 
stricted to the classical sector, and does not lead to QMC 
sign problems. This allows us to employ large-scale QMC 
simulations to study quantum effects on a frustrated spin 
system. In particular, we study the ground state proper- 
ties of the spin- 1/2 easy-axis Heisenberg model on the or- 
thogonal dimer lattice considered by Shastry and Suther- 
land^. Namely, we consider the XXZ Hamiltonian 



H = J E, 



i,3) 



[-A(Sf SJ + S?S?) + SfS*] 



+J' E«u)) [-A(SfSJ + SfS]) + S*S*\ , (1) 

a variant of the model considered in Ref.Qwith ferromag- 
netic transverse exchange (A > 0) and antiferromagnctic 
Ising exchange interactions J, J' > 0. Here, S; denotes 
a spin- 1/2 degree of freedom on site i of the square lat- 
tice, and the first sum extends over all nearest neighbor 
bonds. The second sum runs over a staggered subset of 
the next-nearest neighbor bonds, as indicated in Fig. 1. 

The model in Eq. (1) maintains the frustrated nature of 
the antiferromagnetic Ising interactions, and introduces 
ferromagnetic spin exchange terms. Employing the well 
known mapping between spin-1/2 degrees of freedom and 
hard-core bosons, the model can be mapped onto an ex- 
tended bosonic Hubbard model of hard-core bosons hop- 
ping along the bonds of the Shastry-Sutherland lattice 
and experiencing repulsive interactions proportional to 
the strength of the Ising exchange. Recently, similar 
hard-core boson models have been studied on different 
lattice geometries, and were found to exhibit interest- 
ing order-by-disorder phenomena when quantum fluctu- 
ations lift an extensive ground state degeneracy from the 
Ising limit A = 0, with new quantum phases emerging. 
Examples include a supersolid phase on the triangular 
lattic o 13 ' 14 ! 15 ' 16 , valence-bond- solid o 1 7 ' 18 and a Z2 spin 
liqui d 19 i 20 ' 21 on the kagome lattice, and a U(l) liquid on 
the pyrochlore lattice 2 ^. In the limit of dominating ki- 
netic terms, such models stabilize a superfluid phase on 
both bipartite and non-bipartite lattices. In magnetic 
language, the superfluid corresponds to a transverse fer- 
romagnetic spin alignment, driven by the ferromagnetic 
nature of the transverse spin exchange. For the remain- 
der of the paper, we prefer using the spin language, but 
occasionally find it convenient to also employ the bosonic 
notation. 

As reviewed in the following section, the antiferromag- 
netic Ising model on the Shastry-Sutherland lattice ex- 
hibits regions of extensive ground state degeneracy, sim- 
ilar to the Ising model on the triangular and kagomc 
lattices. Motivated by the above mentioned studies on 
these frustrated geometries, we here assess the effects of 
quantum fluctuations on the classical degenerate ground 
states on the Shastry-Sutherland lattice, and explore the 
phase diagram of the full quantum model. We find in 




FIG. 2: (a) Possible spin configurations in the Ising limit on a 
non-void plaquette for J' > 2 J. (b) Additional configurations 
allowed on a non-void plaquette at J' = 2J. Full (open) 
circles denote spin up (down) states. 



this system a dimer triplet state, discussed in detail be- 
low, to emerge out of the classical degenerate region. In 
addition, the system shows a Neel ordered phase and a 
superfluid regime. We study the quantum phase tran- 
sitions between these different phases, and consider the 
effects of a magnetic field. We do not obtain indications 
for supersolidity in this model, but find that quantum ef- 
fects lead to the stabilization of a magnetization plateau 
at 1/2 of the full saturation, that does not persist down to 
the Ising limit. This is in contrast to the case of e.g. the 
triangular and kagome lattice, where all plateaus found 
in the quantum regime persist down to the Ising limit, 
where they have a largest extension. 

The remainder of the paper is organized as follows: In 
the next section, we review the properties of the antiferro- 
magnetic Ising model on the Shastry-Sutherland lattice. 
Then, we present in Sec. Ill our numerical results on the 
phase diagram of the model introduced above. In order 
to explain in a simple picture the emergence of the dimer 
triplet phase, we employ degenerate perturbation theory 
around the Ising limit, which will be discussed in Sec. 
IV. In Sec. V, wc analyze the properties of the system in 
finite magnetic fields, discuss the appearing magnetiza- 
tion plateaus, and scan for supersolid phases. Finally, we 
conclude in Sec. VI by relating our numerical finding to 
the properties of the isotropic Heisenberg antiferromag- 
net (A = —1) on the Shastry-Sutherland lattice, and 
discuss connections to recent studies on its magnetiza- 
tion process. We also comment on the recent work on 
the compound TmB4, suggested to realize the easy- axis 
antiferromagnetic Shastry-Sutherland model close to the 
Ising limit (-1 < A < 0)2. 



II. 



ISING LIMIT 



Before exploring in detail the phase diagram of the 
quantum spin model introduced above, it is convenient 
to review the properties of the Ising limit, A = 0, dis- 
cussed in Ref. |j. In the Ising limit, the model in Eq. (1) 
stabilizes an antiferromagnetic Neel phase for sufficiently 
weak J', up to J'/J < 2. For J' > 2J the clas- 
sical system has a macroscopically degenerate ground 
state manifold with an extensive ground state entropy 
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S = []n(2)/2]k B N = 0M7k B N : from all configurations 
that cover each of the J' dimer bonds with a pair of oppo- 
site spins. Here, N denotes the number of spins. Exactly 
for J' = 23, the degeneracy of the ground state manifold 
is further enlarged, as additional low-energy configura- 
tions proliferate. Shastry and Sutherland proved a lower 
bound S > 0.4812 k B N on the ground state entropy at 
J' = 2J, via mapping the model to a 10-vertex model 
and using brading techniques^. A simple estimate of the 
ground state degeneracy can be obtained by employing 
the argument from Pauling's estimate of the residual en- 
tropy of ice22. For this purpose, consider one of the filled 
plaquettes on the Shastry-Suthcrland lattice. While for 
J' > 23 the eight configurations shown in Fig. 2 (a) 
provide minimal contributions of this plaquette to the 
total energy, for J' = 23 the two configurations shown 
in Fig. 2 (b) also lead to a minimal energy contribution. 
Given that out of the 16 possible configurations of the 
four spins forming the plaquette these 10 configurations 
are thus feasible, we obtain an estimate for the ground 
state entropy 



S ~ k B In 



-) 

16 J 



N/2' 



0.458 k B N, 



(2) 



to be compared to the above bound by Shastry and 
Sutherland. We note, that for 3' > 23 the Pauling esti- 
mate S = k B ln(2 JV / 2 ) recovers the exact result. 

Besides the Ising limit, Shastry and Sutherland con- 
sidered the effects of antiferromagnetic transverse spin 
exchange terms (A < in our notation), and proved 
that the system possesses an exact dimer-singlet product 
eigenstate, that at least for 3' > 23 becomes the system's 
ground stated. Later studies by various groups consid- 
ered the full quantum phase diagram of this model, which 
up to date is not conclusively established (c.f. Ref. for 
a review of the various theoretical and numerical pro- 
posals), even though numerical evidence has been put 
forward, that the SU(2) symmetric model (A = —1) 
features indeed three phases: (i) a low- J' antifcrromag- 
netically ordered Neel phase, (ii) the large- J' dimer sin- 
glet phase, and (iii) an intermediate valence bond crys- 
tal (VBC) phase, which breaks the lattice symmetry by 
forming resonating plaquette singlet states on one of the 
subsets of the void plaquettes of the Shastry-Sutherland 
lattice 2 ^. This VBC phase has a two-fold degenerate 
ground state and a finite spin excitation gap. For the 
remainder of this work, we study the quantum model in 
the region A > 0, where large-scale QMC simulations 
are possible, in contrast to the previously studied case 
of A < 0. Furthermore, this model relates directly to a 
model of hard-core bosons, as mentioned in Sec. I. 



III. QUANTUM PHASE DIAGRAM 

In this section, we present our numerical results on the 
phase diagram of the model in Eq. (1). These results are 
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FIG. 3: (Color online) Ground state phase diagram of the 
spin- 1/2 XXZ model on the Shastry-Sutherland lattice with 
ferromagnetic transverse spin exchange. The dotted (solid) 
line denotes a first-order (continuous) quantum phase tran- 
sition. Uncertainties on the indicated phase boundaries are 
below the symbol size. The dashed line is a guide to the 
eye indicating a line of constant J'/J = 2. The dashed- 
dotted line gives the estimated phase boundary of the dimer 
triplet phase within perturbation theory around the point 
(J'/J, A) = (2,0), discussed in Sec. IV. 



based on QMC simulations of finite systems with up to 
N = 36 x 36 lattice sites, using period boundary condi- 
tions. In the simulations, we scaled the inverse temper- 
ature as j3 = 1/T = 8L/A3 in order to access ground 
state properties. Here, L denotes the linear system size. 
The QMC simulations were performed employing a gen- 
eralized directed-loop updat e 25 i 26 in the stochastic series 
expansion (SSE) algorithm^!. For the results obtained on 
the larger lattices, in particular in finite magnetic fields, 
we employed a decoupling of the Hamiltonian in plaque- 
tte terms, instead of the more conventional bond decou- 
pling for the SSE formulation. 

In Fig. 3, we present the ground state phase dia- 
gram resulting from our calculations. We find that the 
extend of the antiferromagnetically ordered Neel phase 
shrinks essentially linearly upon increasing A from 
up to the Heisenberg point at (A, 3'/ 3) = (1,0) (for 
3' = 0, the model reduces to a spin model on the bi- 
partite square lattice, and the sign of A can be inverted 
by an unitary transformation , thus relating the point 
(A, 3' j 3) = (1,0) to the isotropic Heisenberg model at 
(A, J'/J) = ( — 1,0)). In hard-core bosonic language, the 
Neel phase corresponds to a checkerboard solid with al- 
ternating occupation of the lattice sites. In our QMC 
simulations, we determine the corresponding structure 
factor Saf for antiferromagnetic order, 



SAF = ^ e%€ 3 ( S i S 



(3) 



h0 



where e, = ±1, depending on the sublattice, to which 
lattice site i belongs. Neel order is present, if in the ther- 
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FIG. 4: (Color online) Spin stiffness ps at fixed A = 1/2 as 
a function of J'/J for different system sizes. 



modynamic limit Saf/N scales to a finite value. For 
dominant transverse exchange, A > 1, the model re- 
duces to an ferromagnetic XY model on the Shastry- 
Sutherland lattice, which in bosonic language relates to 
an non-frustrated tight-binding hopping model. Hence, 
the system is expected to exhibit a bosonic superfluid 
phase for large values of A > 0, which in spin language 
relates to a ferromagnetic ordering within the XY plane. 
Such a phase is characterized by a finite value of the su- 
perfluid density, or spin stiffness, which in the QMC sim- 
ulations can be obtained from measuring the spin wind- 
ing number fluctuation^ (W 2 ) as 



Ps 



T 
AJ 



(W 2 ). 



(4) 



As an example, we show in Fig. 2] the behavior of ps 
for different system sizes as a function of J'/J for fixed 
A = 1/2. A region with finite spin stiffness is found 
for 0.9 < J'/J < 2.5. The strong discontinuity of ps 
at J'/J = 0.9 indicates that the quantum phase transi- 
tion between the Neel ordered phase and the superfluid 
is strongly first order. Such behavior is also seen by mon- 
itoring the antifcrromagnctic structure factor Saf upon 
crossing the phase boundary, as shown in Fig. [5l again 
for A = 1/2. Combining the results for Saf and ps, we 
obtain no indication for an intermediate region exhibit- 
ing both finite superfluidity and diagonal long-range or- 
der as inside a supersolid phase, as expected from the 
commensurate half-filling of the lattice. For dominant J' 
(e.g., J'/J > 2.5 at A = 1/2), both Saf and ps eventu- 
ally vanish. We explicitly verified that inside this regime 
the model does not exhibit long-ranged correlations in 
the longitudinal nor the transverse spin-spin correlation 
function. In addition, also the bond-order-wave struc- 
ture factors do not exhibit long-ranged order in the spin 
exchange correlation function (corresponding to kinetic 
energy correlations in the bosonic model). 

Since the parameter region with J' > 2J approaches 
the degenerate region of the Ising limit for A — > 0, quan- 
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FIG. 5: (Color online) Antiferromagnetic structure factor 
across the phase boundary between the Neel ordered phase 
and the superfluid regime at fixed A = 1/2 as a function of 
J'/J for different system sizes. 



turn effects indeed select a unique phase from this de- 
generate ground state manifold. In particular, for small 
values of A <C 1, the ground state in this large- J' regime 
can be obtained using degenerate perturbation theory 
around the Ising limit, discussed in detail in the next 
section. Within first-order perturbation theory in A one 
then finds that for J' > 2 J quantum fluctuations select 
the following dimcrized state of localized S% ot = triplet 
states on each dimer: 



|^) = ®^(IU> + Ut» d . 



(5) 



Here, the direct product extends over all J' dimer bonds 
on the lattice. Obviously, this symmetric linear local 
combination results from the ferromagnetic nature of the 
transverse spin exchange (A > 0) considered here. For 
A < 0, one instead recovered the exact dimer singlet 
state found by Shastry and Sutherland^. In contrast to 
the dimer singlet state however, the above state \iPd) is 
not an eigenstate of the Hamiltonian for finite values of J. 
As discussed in the following sections, processes in higher- 
order perturbation theory lead to local correlations be- 
tween the dominant resonances on the dimers. Hence, 
the ground state in the large- J' region of the quantum 
phase diagram does not take the above direct product 
form, but approaches it for J/ J' — > 0. From the ground 
state energy, we still find that the state \iPd) provides a 
appropriate variational state for the true ground state in 
the large- J' regime. This can be seen even at A = 1/2, 
i.e. significantly away from the Ising limit, from a com- 
parison between the system's ground state energy E and 
the variational energy of \iPd), 

£ D HMffhM = -^p(i + 2A), (6) 

shown in Fig. [SJ Due to the dominant formation of 
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FIG. 6: (Color online) Comparison of the ground state energy 
E as a function of J'/ J for A = 1/2 and the variational energy 
of the dimerized state of localized S^ ot = triplet states, 
E D = (ipn\H\i> D ). 

Sf ot = triplet states on the dimer bonds, we denote this 
magnetically disordered phase as a dimer triplet phase. 

Since no spatial symmetry is broken in the dimer 
triplet phase, we expect the quantum phase transition 
from the superfluid with broken U{1) symmetry to the 
dimer triplet phase to be continuous, and to belong in 
the universality class of the three-dimensional (3D) 0(2) 
model, with a dynamical critical exponent z = 1. In 
order to study the nature of this quantum phase transi- 
tion in the QMC simulations, we scanned the transition 
region at fixed values of either A or J' /J, varying the 
other parameter through the phase boundary. Denot- 
ing the varied parameter by X, at a continuous quantum 
phase transition the spin stiffness scales as 

p s (X,L)=L-*f(t x L 1 /»,P/L*) (7) 

with a scaling function /, and the correlations length 
exponent v. Here, 

X — X r , . 

tx = — ^ (8) 

denotes the relative distance away from the critical point 
at X = X c . From the above scaling relation, it follows 
that X c can be determined as the crossing point of finite 
size data for the rescaled spin stiffness L z ps- Further- 
more, with appropriate values of the critical exponents 
z and ^, the scaling function f(-,A) is then obtained by 
plotting L z ps vs. txE x l v for a fixed value of f3/L z = A. 
As an example, we consider a scan in X = J' / J at a fixed 
value of A = 1/2, for which the finite size data of ps is 
shown in Fig [7] (a). For z = 1, we obtain X c = 2.46(2) 
from a clear crossing point in Fig [7] (b) , and a clear data 
collapse within a finite critical region, taking v = 0.6723 
for the 3D 0(2) models, as shown in Fig. [7J (c). 

Proceeding this way for other values of A, we eventu- 
ally obtained the phase boundary shown in Fig. 3. From 




FIG. 7: (Color online) Spin stiffness ps at fixed A = 1/2 as 
a function of J' /J for different system sizes (a), and rescaled 
with linear system size L (b), with J' /J — 2.46 marked by 
the dashed line. Part (c) shows the data collapse expected 
from a finite size scaling analysis for the 3D 0(2) universality 
class. 



this analysis, we find that at low values of A < 0.225, the 
dimer triplet phase extends below the line J' / J = 2, in- 
dicated by the dashed line in Fig. 3. In order to illustrate 
this explicitly, we show in the left panels of Fig. [8] the re- 
sults of the finite size scaling analysis at A = 1/6, where 
the transition point is located at J 1 / J = 1.928(5) < 2. 
As a crosscheck, the right panel of Fig. [5] shows the data 
of the finite size scaling analysis at fixed J'/ J = 1.928, 
where the transition is indeed observed at A = 1/6. Sim- 
ilarly, when varying A at fixed J' /J = 2 a transition 
point between the dimer triplet phase and the superfluid 
regime is found at A = 0.225(1), as extracted from Fig. [9] 
These results suggest that either (i) the superfluid region 
separating the Neel ordered phase and the dimer triplet 
phase persists down to the Ising limit, or (ii) the first 
order transition line and the second order transition line 
meet at a finite value of A, or (iii) an additional phase 
appears near J' = 2J for even smaller values of A < 0.1. 

Such an additional phase might be expected to be se- 
lected by quantum effects from the state of enlarged de- 
generacy of the Ising model at J' = 2J for finite A. 
The QMC simulations could not be extended to signifi- 
cantly smaller values of A, due to an reduced efficiency in 
parameter regions dominated by the frustrated diagonal 
part of the Hamiltonian. However, as discussed in the 
following section, degenerate perturbation theory in A 
indicates that at J' = 2 J the dominant effect of a finite 
A is to effectively drive the system away from J' = 2 J 
towards the region J' > 2 J. This leads us to exclude op- 
tion (iii) from the above list. Scenario (ii) would imply 
a direct first order transition between the Neel ordered 
phase and the dimer triplet state for sufficiently low val- 
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FIG. 8: (Color online) Left panel: Spin stiffness ps at fixed 
A = 1/6 as a function of J'/J for different system sizes (a), 
and rescaled with linear system size L (b), with J'/J = 1.928 
marked by the dashed line. Part (c) shows the data collapse 
for the 3D 0(2) universality class. Right panel: Spin stiffness 
ps at fixed J'/J = 1.928 as a function of A for different 
system sizes (d), and rescaled with linear system size L (e), 
with A = 1/6 marked by the dashed line. Part (f) shows the 
data collapse for the 3D 0(2) universality class. 



ues of A > 0, whereas within scenario (i) the superfluid 
phase would always separate the two phases. 



IV. PERTURBATION THEORY 

In order to study more closely the emergence of the 
dimer triplet phase from the Ising limit upon introduc- 
ing transverse exchange interactions, we employed degen- 
erate perturbation theory, starting from the degenerate 
ground state manifold in the Ising limit A = 0. First, 
we consider the region J' > 2 J, where the degenerate 
ground state manifold is spanned by independently plac- 
ing two opposite spins on each J'-dimer, as discussed in 
Sec. II. 

For each such J'-dimer d, we denote these two lowest 
energy states as 



a)d = I u>d, 

ii)d = I IT)d, 



(9) 



which form an effective spin- 1/2 degree of freedom on the 
dimer d. We separate the full Hilbert space of the system 
into the model space M , spanned by these ground state 
configurations, and the orthogonal space O. A basis of 
M is given by the orthonormal states 



\^ a ) d , with|^) rf G{|^) d ,|ir) d }- (io) 



The orthogonal space O is spanned by all states of the 
Ising model, that do not belong to this set. We denote 
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FIG. 9: (Color online) Spin stiffness ps at fixed J'/J = 2 
as a function of A for different system sizes (a), and rescaled 
with linear system size L (b), with A = 0.225 marked by the 
dashed line. Part (c) shows the data collapse expected from a 
finite size scaling analysis for the 3D 0(2) universality class. 



these orthonormal basis states by \ipb), for which at least 
one dimer d has both spins equal, i.e. | ]])d or | ll)<2- 
The Hamiltonian of Eq. (1) is similarly split into the 
model Hamiltonian Hq and a perturbation part H± oc A, 



H = 



(id) 



J' 



{(id)) 



SfS* 



J y £[-A(3fS^ + 3¥S«)] 

(id) 

((id)) 



(11) 



where Hq is diagonal in the basis of M and O introduced 
above. The effective Hamiltonian H e ff, that describes 
the effective dynamics induced by Hi within the model 
space M is given by degenerate perturbation theory up 
to third order in A as^ 

H efJ = PH Q P + PHiP^ + PH1RH1P 

1st order 2nd order 

+ PH1RH1RH1P - PHi RRHi PHiP 



3rd order 



+0(A 4 ), 



where 



(12) 



(13) 



is the projection operator onto the model space in terms 
of the above constructed basis states \ip a ), and 
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FIG. 10: (Color online) Spin exchange processes contributing 
to the effective Hamiltonian H^jj in first order perturbation 
theory. Full (open) circles denote spin up (down) states. 



is the resolvent operator with Eq = —NJ'/8 the degen- 
erate ground state energy of H , and Eq = (V , &|i?o|V'&) 
the energy of the basis state \ifjt>) of the orthogonal space 
O. We start by considering the first order contribution 
to H e f f . The only possible process in this case is a single 
spin-flip along a diagonal bond, graphically represented 
in Fig.[l0l In terms of effective spin operators and 

Sji, that act on the effective dimer spin states | | JJ-)d, 
the effective Hamiltonian in first order perturbation thus 
reads, 



H 



(i) 

eff 



(15) 



corresponding to a uniform transverse magnetic field act- 
ing on the effective dimer spins. The lowest-energy eigen- 
state of H e f f for A > is the direct product state 



1 

71 



(\itu+\m 



i 



V2 



(IU) + UT)) d , 



referred to already in Eq. |(SJ) of Sec. Ill, corresponding 
to the decoupled dimer state with each dimer forming a 



Of, 



triplet state. In case of an antiferromagnetic 



transverse exchange, A < 0, the lowest energy state of 
H e ff is the dimer singlet state 

= ® -4 (i tu - 1 m = 4 (i u> - 1 it» d , 

d V d V 

proven to be an exact eigenstate of the full Hamiltonian 
H for A < by Shastry and Sutherland^. However, \i))d) 
is not an eigenstate of the full Hamiltonian, and correla- 
tions between the effective dimer spins are introduced in 
higher order perturbation theory. 

In order to assess the nature of these correlations, we 
consider processes occurring in second order perturbation 
theory. These involve two spin exchanges along axial 
bonds. Apart from diagonal terms, the result of such 
processes is again to flip the effective spin on one of the 
dimers, as shown in Fig. [TTJ The matrix element of each 
such process depends in detail on the specific local spin 
configuration on the dimers neighboring the dimer that 
undergoes the spin flip. Among the various possibilities, 
one with the largest amplitude is shown in Fig. [11] (a). In 





FIG. 11: (Color online) Two different spin exchange processes 
contributing to the effective Hamiltonian in second or- 

der perturbation theory. These processes mediate the forma- 
tion of an Sf 0i = triplet state on the corresponding dimer, 
while blocking its formation on the neighboring dimers (inside 
ellipses). Full (open) circles denote spin up (down) states. 
Numbers indicate the order of the spin exchange in the pro- 
cesses. 



(2) 

the second order effective Hamiltonian H f f , this process 



contributes a term 



s z k ){\ 



Sf 



(16) 



which provides a further contribution to the transverse 
field operator at site d, dressed by diagonal operators 
from the neighboring dimers, that project out the specific 
configuration of dimer spin states according to the con- 
figuration shown in Fig. QT] While the transverse field in 
H^X acts locally on each dimer, the dressed transverse 

(2) 

field operators deriving from H e f $ lead to correlations 
among nearest neighbor dimer spins. For example, the 
spin exchange process on dimer d shown in Fig. [11] (a) 
could not take place as indicated by the arrows, if dimer 
k was in the opposite spin state. Instead, a different pro- 
cess could take place, as shown in Fig. [IT] (b), that leads 

(2) 

to a similar term in H e jp but with a different energy 
denominator than in Eq. (| 16[) . In this way, details of the 
local dimer configuration enter the effective Hamiltonian 
in a rather complex manner. 

The explicit form of the total effective Hamiltonian up 
to second order in A involves several terms, containing 
products of up to five effective spin operators, such as 
the term given explicitly in Eq. ([Hi]) . While the ground 
state of this effective Hamiltonian is not directly accessi- 
ble, the general structure of these terms indicate that it 
will be a dressed version of \iPd)> with local inter-dimcr 
correlations induced by the above virtual spin exchange 
processes. This reflects the QMC result, that the true 
ground state is close in energy to {i/jd}, and does not 
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FIG. 12: (Color online) Spin exchange process contributing 
to the effective Hamiltonian in third order perturbation 

theory. This process leads to a flip of the effective spins on two 
neighboring dimers (inside boxes). As a result of this process, 
the effective spins inside the boxes have been exchanged. Full 
(open) circles denote spin up (down) states. Numbers indicate 
the order of the spin exchange processes. 



exhibit long ranged correlations. 

The effective Hamiltonian up to second order in pertur- 
bation theory does not include transverse spin exchange 
terms, that would contain products of off-diagonal effec- 



tive spin operators such as S~J S k 



S d Sfc on two neigh- 



boring dimers d and k. However, such terms appear in 
third order of the degenerate perturbation theory, e.g. 
from the process shown in Fig. [T^l Similar to the trans- 
verse magnetic field operator in second order, the matrix 
elements depend on the details of the spin configuration 
on the neighbors of the two dimers that undergo an ef- 
fective spin exchange. The transverse effective spin ex- 
change operators are thus similarly dressed by additional 
projection operators. Becoming more relevant at large 
values of A, we expect that such exchange terms even- 
tually drive the transition from the dimer triplet phase 
to the superfluid region, that was observed in the QMC 
simulations (c.f. Fig. [3]). 

We now turn to the special point J' ~ 2 J, where the 
ground state in the Ising limit has an enhanced degener- 
acy, as discussed in Sec. II. Consider one of the J'-dimcr 
of the Shastry-Sutherland lattice in the Ising limit A = 0. 
For J' = 2 J, this dimer can be either in one of the local 
configurations of Fig. 2 (a), which are also allowed con- 
figurations for J' > 2J, or it is part of one of the local 
configurations shown in Fig. 2 (b). Introducing a finite 
A > in first order perturbation theory, a splitting in 
the energy between the configurations of Fig. 2 (a) and 
those of Fig. 2 (b) results, since the configurations of 
Fig. 2 (b) cannot gain exchange energy by a transverse 
spin exchange along the dominant J' bonds, in contrast 
to the configurations of Fig. 2 (a). To first order in A, 
this energy difference equals 



SE = -J + J - + ^. 

2 2' 



(17) 



near the point (A, J' /J) = (0,2). A finite SE > thus 
leads to a partial lifting of the ground state degeneracy, 
and only the local configurations of Fig. 2 (a) remain to 
span the low-energy sector. The configurations of Fig. 2 



(b) are split-off by a finite energy difference SE > from 
the ground state manifold. Since this energy difference 
remains positive for 



SE > 



,r_ 2 

J > 1 + . 



(18) 



we expect that as long as J' / J > 2/(1 + A) close to 
(A, J'/ J) = (0,2), the system is driven by quantum ef- 
fects towards the same phase as for J' > 2J. In fact, 
this expectation is in agreement with the QMC phase di- 
agram of Fig. 3, where we found that (i) for small A > 
at J' = 2J, the system enters the dimer triplet state 
as it does for J' > 2 J, and (ii) the phase boundary of 
the dimer triplet phase closely follows the limiting line 
J'/ J = 2/(1 + A) according to SE = for small A (com- 
pare to the dashed-dotted line in Fig. 3). The above 
argument was based on energy considerations on an iso- 
lated dimer. Due to this restriction, we are not able here 
to discern, if the superfluid phase indeed terminates at 
a small, but finite values of A, or if it persists down to 
any finite value of A > 0. This would require the intra- 
dimer exchanges to be taken into account within higher 
orders of perturbation theory. However, from our analy- 
sis of the case J' > 2J discussed above, we expect that 
also in this case, the effective model will not allow for an 
explicit solution, thus leaving this question unanswered. 
Hence, here we do not attempt to extend on this issue, 
but instead move on to study the magnetization process 
in the model under consideration. 



V. MAGNETIZATION PROCESS 

After exploring the ground state phase diagram, we 
now consider the model of Eq. (1) in the presence of a 
finite magnetic field h, which couples to the spins by the 
standard Zeeman term, 



H 



(19) 



From considering a singly flipped spin with respect to 
the fully polarized state, one finds that the system is 
fully polarized for magnetic fields h beyond 



h s = (2J + J'/2)(l + A). 



(20) 



Before discussing the magnetization process of the full 
quantum model, it is again useful to consider first the 
Ising limit. The authors of Ref. state that the Ising 
model on the Shastry-Sutherland lattice exhibits a mag- 
netization process with a single plateau at m/m s = 1/2, 
where m denotes the magnetization, and m s its satura- 
tion value, at least for J' < 2J. In our notation, this 
1/2 plateau should extend between 2 J — J'/2 < h < 
2J + J'/2. In particular, in case J' = J, well within 
the Neel ordered zero-field regime, this range becomes 
3/2 < h/J < 5/2. This conclusion in Ref. [9| was based 
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FIG. 13: (Color online) Magnetization curves m/m s as func- 
tions of h for the Ising model on the Shastry-Sutherland lat- 
tice at J' = J for different system sizes, obtained from Monte 
Carlo simulations. The inset shows the results for the 4x4 
and 8x8 lattices. 



FIG. 14: (Color online) Magnetization curves m/m s as func- 
tions of h for the Ising model on the Shastry-Sutherland lat- 
tice at J' = 2.47 for different system sizes, obtained from 
Monte Carlo simulations. 



on analyzing a finite system with 16 spins only. In or- 
der to check this conclusion on larger system sizes, we 
performed a systematic finite size analysis, using classi- 
cal Monte Carlo (MC) simulations on lattices with up 
to 18 x 18 spins in the canonical ensemble. For these 
simulations, we employed a single spin- flip Metropolis 
algorithm, and allowed for a simulated annealing of the 
system from a large initial temperature T ~ J down to 
the final temperature during an initial state of equilibra- 
tion. This way we were able to obtain MC results down 
to T/J = 0.1 for L = 4,6, and T/J = 0.4 for L = 18. 
While for the purpose of the current study, one can draw 
relevant conclusions on the magnetization process from 
these data, it will be interesting to obtain more refined 
numerical data on the magnetization process in the Ising 
limit, using extended ensemble sampling methods, simi- 
lar to the approach taken in Ref. |3jJ for the square and 
triangular lattice. However, this lies beyond the scope of 
the current study, which is directed towards the quantum 
regime. 

Upon presenting the MC results, we first discuss the 
case J' = J considered in Ref. 0. Our MC data for the 
magnetization process on different lattices are collected 
in Fig. [T3l The inset of Fig. [13] shows our data for a 
4x4 system, which appear to confirm the conclusion of 
Ref. i. However, upon increasing the system size, we 
find that different plateau structures appear. In general, 
in order to allow a system to establish a certain magne- 
tization plateau, appropriate lattice sizes and boundary 
conditions must to be chosen. Otherwise, geometric con- 
straints could frustrate certain magnetization patterns. 
In the present case, it can be seen from the numerical 
data, that only for L a multiple of 3, the system estab- 
lishes a wide m/m s = 1/3 plateau, which is consistently 
observed for L = 6,12 and 18. Instead, the L = 4 and 
8 systems cannot establish the corresponding magnetic 
superstructure, and hence lead to rather different magne- 



tization curves with strong finite size effects, that do not 
represent thermodynamic limit behavior. The magneti- 
zation curve in Fig.[T3]for L = 6, taken at T/J = 0.1, ex- 
hibits a magnetization plateau at m/m s = 1/3, extend- 
ing from h = J up to the saturation field at h = 5/2 J. 
The data shown for the two larger systems, for which 
higher temperatures had to be taken in the MC simu- 
lations, are consistent with a thermal smothering of the 
magnetization jumps out of the plateau towards to = 
and full saturation, respectively. We conclude that for 
,/' = J the Ising model exhibits a single intermediate 
magnetization plateau at m/m s = 1/3, extending from 
h = J up to the saturation field at ft. = 5/27. The 
different conclusion of Ref. appears to be due to the 
usage of inappropriate finite lattice sizes. Note, that all 
values of L considered here were even, i.e. a plateau at 
m/m s = 1/2, if it would exist in this model, would not 
be frustrated by finite lattice effects. In fact, as discussed 
below, we find that such a 1/2 plateau appears for finite 
values of A due to quantum effects. Next, we consider 
the case J' = 2.4J, well inside the degenerate region of 
J' > 2J. Below, we will compare the Ising model re- 
sult to QMC data on the magnetization process for finite 
A > at the same value of J'/ J. The MC results for 
the magnetization process in the Ising limit are shown 
in Fig. 1141 where we now consider linear system sizes 
L that arc a multiple of 6. Again, we observe a wide 
1/3 magnetization plateau, which extends from h m 0.87 
up to the saturation field at h = h s = 3.2J. A precise 
estimate of the lower boundary of the plateau is not ac- 
cessible from the current finite size data, as we could not 
collect data at sufficiently low temperatures on larger sys- 
tems. There appears however a finite magnetization with 
a smooth increase well before the plateau is entered. It 
will be interesting to explore this low— m regime in more 
detail using extended ensemble methods. The point that 
is important for the following discussion, and which fol- 
lows also from the current MC data, is the absence again 



10 




h/J 

FIG. 15: (Color online) Magnetization rn/m 3 and spin stiff- 
ness ps for J' — 2. 4 J and A = 1/4 as functions of the applied 
magnetic field h. 



of a magnetization plateau at m/m s = 1/2. Instead, the 
magnetization exhibits a jump from the 1/3 plateau up 
to magnetic saturation at ft. = h s . 

As we show next, a plateau at m/m s = 1/2, while 
not obtained in the Ising limit, emerges in the quan- 
tum model at finite values of A > 0. Fig. [15] shows 
the magnetization process for J' = 2. 4 J and A = 1/4 
from QMC simulations. We again find a magnetiza- 
tion plateau at m/m s = 1/3 as well as a plateau at 
m/m s = 1/2. Adding a transverse exchange to the Ising 
model thus leads to a softening of the large magnetiza- 
tion jump from m/m s = 1/3 to 1 in the Ising limit, and 
an additional plateau region appears with m/m s = 1/2. 
In hard-core bosonic language, the transverse exchange 
maps onto a non-frustrating hopping amplitude. In the 
current situation, this finite boson hopping does not lead 
always to superfluidity, but drives the system into an 
insulating phase at filling p = 3/4, corresponding to 
m/m s = 1/2 (due to particle- hole symmetry, a similar 
insulating region emerges also for p = 1/4, correspond- 
ing to m/m s = —1/2). Indeed, we find from Fig. 1151 
that the spin stiffness p s vanishes inside both plateau 
regions. It is then in order to study, if the magnetic 
excitations in both plateau phases form period crystals, 
and what the structures of such solid arrays would be. 
Long-range crystalline ordering of magnetic excitations 
inside magnetic plateau regions has previously been an- 
alyzed for the isotropic spin- 1/2 Heisenberg model on 
the Shastry-Sutherland model^^^^^^^^, based 
on different approximative schemes. The relevant struc- 
tures expected from these studies are shown in the inset 
of Fig. for the 1/3 (upper panel), and the 1/2 (lower 
panel) plateau, respectively. They are expressed in terms 
of periodic arrangements of S^ ot = 1 dimcr triplet states 
denoted by dumbbells, in a background of less po- 
larized, e.g. S£ ot = states on the remaining dimcrs. We 
now assess, if these crystalline patterns appear also in the 
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FIG. 16: (Color online) Dimer super structure factors 
5(2^/3,2^/3) (upper panel) and S(n,n) (lower panel) at 
J' — 2. 4 J and A = 1/4 as functions of the applied mag- 
netic field h. The insets illustrate the corresponding magnetic 
superstructures, where a dumbbell denotes a fully polarized 
dimer. Nd = N/2 is the number of dimers. 



current quantum model, which represents an anisotropic 
version of the Heisenberg model considered thus far. 

In a QMC simulation, one can probe exactly for these 
magnetic superstructures by measuring an appropriate 
structure factor. Hence, we analyze the ordering pat- 
tern of the magnetic excitations by measuring the triplet 
excitation structure factor 



5 ^) = ]v-E ei;q(rd " rfc) ^)' ( 21 ) 

d d.k 



where Pd is a projector on the S^ ot = 1 dimer triplet 
state | |t)d on dimer d. Nd = N/2 equals the number 
of dimers in the finite system of N sites. The positions 
rd of the dimers and the momentum space vector q are 
defined with respect to the convenient coordinate system 
formed by the square lattice of dimers, with the distance 
between the centers of two neighboring dimers taken as 
unity. In this way, the solid order shown in the upper 
panel of Fig. \W\ for m/m s = 1/3 corresponds to a peak 
in S(q) oc Nd at q = (27r/3, 27r/3), and at q = (n, n) 
for the order shown in the lower panel of Fig. [1^] for 
m/m s = 1/2. The main panels of Fig. fTBl show these 
components of the structure factor for different system 
sizes, respectively. We checked, that no other signal in 
5(q) appeared, apart from those explicitly shown. From 
these data we conclude, that indeed the crystalline or- 
ders of the magnetic excitations shown in Fig. \W\ are 
stabilized within the two magnetization plateau regions, 
respectively. In order to study the real-space distribu- 
tion of the magnetization among the lattice sites in more 
detail, we also measured the mean local values of the 
magnetization. The obtained distributions of the local 
magnetization are shown for the 1/3 and 1/2 plateau in 
Fig. [13 (a) and (b), respectively. They agree well with the 
distributions reported for the isotropic Heisenberg model 
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FIG. 17: (Color online) Distribution of the local magnetiza- 
tion mi/iris for the magnetization plateaus at m/m s = 1/3 
(a) and 1/2 (b), represented as indicated by the size of circles 
on the lattice sites. Results are shown for a 24 x 24 lattice. 



in the studies mentioned above. 

Upon varying h beyond these wide magnetization 
plateaus, the magnetization undergoes discontinuous 
jumps, as seen in Fig. 1151 alert when entering the 1/2 
plateau from below, where such a jump appears not that 
well resolved. The spin stiffness p s exhibits a similar 
behavior, with clear jumps upon entering the plateaus, 
expect when entering the 1 /2 plateau from below. Apart 
from this case, the supcrfiuid-insulator transitions are 
thus clearly first order. Concerning the region below the 
1/2 plateau, one might instead consider the presence of 
a supersolid state of magnetic excitations, with both a 
finite superfluid density and a crystalline superstructure. 
In order to assess, if supersolid behavior is indeed present 
in this regime, we show in Fig.[T5]thc finite size scaling of 
both the superfluid density p s and the relevant structure 
factors from the neighboring magnetization plateaus at a 
magnetic field of h = 2.6, inside the possibly supersolid 
region. From the finite size scaling we can exclude the 
presence of a supersolid phase in this parameter regime. 
Between the two magnetization plateaus, the system is 
thus in a uniform superfluid state. Our analysis of the 
magnetization process at J'/ J = 2.4 and A = 1/4 did 



FIG. 18: (Color online) Finite size scaling of the superfluid 
density ps and the dimer superstructure factors of the neigh- 
boring plateaus at h = 2.6 for J' — 2. 4 J and A = 1/4. 



not exhibit the presence of any additional plateaus at 
lower values of m. Indeed, at low-m the magnetization 
curve as well as the superfluid density appear smooth 
in Fig. 1151 In particular, we did not find any of the 
fractional magnetization plateaus mentioned in Sec. I 
for the compounds SrCuB 2 (B0 3 ) 2 and TmB 4 , nor those 
found for the isotropic spin-1/2 Heisenberg model on the 
Shastry-Suthcrland model in recent studies ^ 32 ' 33 . We 
postpone a discuss on this issue to the concluding sec- 
tion. 

Finally, we present our results for the magnetization 
process upon varying over a wider range of parameters 
space. The results of QMC simulations similar to those 
discussed in detail above are summarized in Fig. 1191 
which shows the magnetic phase diagram for a generic 
fixed ratio of J'/ J = 2.4. The figure shows both the 
m/m s = 1/3 magnetization lobe that extends from the 
Ising limit up to A = 0.32(2), as well as the emergent 
m/m s = 1/2 lobe, that extends up to a similar value of 
A, but shrinks upon approaching the Ising limit, where 
no 1/2 plateau persists. 



VI. CONCLUSIONS 

We studied the ground state phase diagram and the 
magnetization process of the spin-1/2 easy-axis XXZ 
model of Eq. (1) with ferromagnetic transverse spin ex- 
change on the Shastry-Suthcrland lattice. The model 
exhibits a Neel ordered phase for small values of J' / J 
and A, and a superfluid phase for dominant A. For suf- 
ficiently large values of J' / J, a dimerized phase with the 
dominant formation of Sf ot = triplet states on the J'- 
dimer bonds is stabilized, that connects to the degenerate 
phase of the model in the Ising limit. We like to compare 
this findings to the case of the isotropic Heisenberg model 
on the same lattice. Also in this model there is strong 
evidence for an intermediate phase 3 -, separating the Neel 
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FIG. 19: (Color online) Magnetic phase diagram of the spin- 
1/2 XXZ model on the Shastry-Sutherland lattice with ferro- 
magnetic transverse spin exchange at J'/ J = 2.4 depending 
on the anisotropy A and the applied magnetic field h. 

ordered phase from the dimer singlet phase. In contrast, 
this intermediate phase is however not a supcrfluid, but 
appears to be realize a valence bond crystal driven by the 
frustrated nature of the transverse exchange interactions 
in this regional. 

Using classical Monte Carlo simulations to assess the 
magnetization process in the Ising limit, we found in con- 
trast to previous claims^ no magnetization plateau at 1 /2 
of full saturation, but instead a 1/3 plateau, with a jump 
towards full saturation. The quantum model shows the 
presence of a 1/2 plateau in addition to the 1/3 plateau. 
This emergence of a magnetization plateau upon adding 
quantum fluctuations (via a finite A) to the Ising model 
remains of the previously observed emergence of a mag- 
netization plateau via thermal fluctuations in the classi- 
cal Heisenberg model on the kagome lattice^. While we 
examined the case of a ferromagnetic transverse spin ex- 
change, the experimental observation of a 1/2 plateau in 
the compounds TmB4 and SrCuB2(B03)2 can be taken 
as indication, that such a plateau also emerges for antifer- 



romagnctic transverse exchange. Indeed, a 1/2 plateau 
is consistently reported also in the studies on the mag- 
netization process of the isotropic Heisenberg model on 
the Shastry-Sutherland lattice 3 -. 

With respect to other reported magnetization plateaus 
both from experiments and in recent theoretical work, we 
did not obtain evidence from our quantum Monte Carlo 
simulations, that they are stabilized in the current model. 
In particular, the pronounced 1/6 plateau, consistently 
reported in recent theoretical wor k 5 ' 32 ! 33 , has not been 
found for the parameter region we examined in our sim- 
ulations (we performed various scans inside the range 
2 < J' I J < 3, down to A = 1/7, with results similar 
to those presented in detail above), even though the pro- 
posed magnetic unit cell 5 - ' 32 ' 33 is commensurate with the 
finite lattices that we employed in our simulations. We 
explicitly checked that no signal in the relevant compo- 
nent of S(q) at q = (7r,7r/3) appears near m/m s = 1/6 
in this model. Hence, we conclude that the model consid- 
ered here exhibits magnetization plateaus at 1/2 and 1/3 
of the full saturation, only. Future work could explore in 
more detail the magnetization process in the Ising model 
inside the low-field region using extended ensemble meth- 
ods, similar to the approach taken in Ref. |3l|. It would be 
interesting to apply the theoretical approaches employed 
in Refs. [32j and [33| to the current model, in order to to 
check their predictions against the unbiased large-scale 
quantum Monte Carlo results present here. 
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